/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include "lower_bound.h"
#include "symmetry.h"
#include "protocol.h"
#include "compatibility.h"

#include "MurmurHash3.h"

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdatomic.h>

typedef struct precomputed_table {
  size_t signature_size;
  size_t num_quads;
  size_t num_hexes;

  size_t signatures_num_bytes;
  uint8_t *signatures;

  size_t meshes_num_bytes;
  uint8_t *meshes;

  size_t minima_num_bytes;
  uint8_t *minima;

  size_t hexahedra_num_bytes;
  uint8_t *hexahedra;

  bool loaded;
} precomputed_table;

static precomputed_table lower_bounds = {.loaded = false};
static precomputed_table known_solutions = {.loaded = false};

static
error_code precomputed_table_load(precomputed_table *table, const char *filename);

static
void precomputed_table_release(precomputed_table *table);

error_code load_lower_bounds(const char *filename) {
  return precomputed_table_load(&lower_bounds, filename);
}

void release_lower_bounds(void) {
  precomputed_table_release(&lower_bounds);
}

error_code load_known_solutions(const char *filename) {
  return precomputed_table_load(&known_solutions, filename);
}

void release_known_solutions(void) {
  precomputed_table_release(&known_solutions);
}

static
error_code load_array(FILE *io, uint8_t **output, size_t *ouput_size);

static
error_code precomputed_table_load(precomputed_table *table, const char *filename) {
  FILE *io = fopen(filename, "r");
  if (!io) return IO_ERROR;

  error_code status = SUCCESS;

  uint8_t buffer[sizeof(uint64_t)];

  if (fread(buffer, sizeof(buffer), 1, io) != 1) {
    status = IO_ERROR;
    goto fail_read_header;
  }

  table->signature_size = decode_uint64(buffer);

  if (fread(buffer, sizeof(buffer), 1, io) != 1) {
    status = IO_ERROR;
    goto fail_read_header;
  }

  table->num_quads = decode_uint64(buffer);

  if (fread(buffer, sizeof(buffer), 1, io) != 1) {
    status = IO_ERROR;
    goto fail_read_header;
  }

  table->num_hexes = decode_uint64(buffer);

  if ((status = load_array(io, &table->signatures,
                           &table->signatures_num_bytes)) != SUCCESS)
    goto fail_load_signatures;

  if ((status = load_array(io, &table->meshes,
                           &table->meshes_num_bytes)) != SUCCESS)
    goto fail_load_meshes;

  if ((status = load_array(io, &table->minima,
                           &table->minima_num_bytes)) != SUCCESS)
    goto fail_load_minima;

  if (table->num_hexes != 0) {
    if ((status = load_array(io, &table->hexahedra,
                             &table->hexahedra_num_bytes)) != SUCCESS)
      goto fail_load_hexahedra;
  }
  else
    table->hexahedra = NULL;

  table->loaded = true;

  if (status != SUCCESS) free(table->hexahedra);
fail_load_hexahedra:
  if (status != SUCCESS) free(table->minima);
fail_load_minima:
  if (status != SUCCESS) free(table->meshes);
fail_load_meshes:
  if (status != SUCCESS) free(table->signatures);
fail_load_signatures:

fail_read_header:
  fclose(io);
  return status;
}

static
error_code load_array(FILE *io, uint8_t **output, size_t *output_size) {
  uint8_t buffer[sizeof(uint64_t)];
  if (fread(buffer, sizeof(buffer), 1, io) != 1)
    return IO_ERROR;

  *output_size = decode_uint64(buffer);
  *output = malloc(*output_size);
  if (!*output)
    return OUT_OF_MEMORY;

  if (fread(*output, *output_size, 1, io) != 1) {
    free(*output);
    return IO_ERROR;
  }

  return SUCCESS;
}

static
void precomputed_table_release(precomputed_table *table) {
  if (!table->loaded)
    return;

  table->loaded = false;

  free(table->hexahedra);
  free(table->minima);
  free(table->meshes);
  free(table->signatures);
}

static
bool compare_boundary(
  solver *solver, quad_index *quads,
  const vertex_index *ref, uint32_t num_quads,
  const uint8_t *valence,
  vertex_index *permutation, vertex_index *inv_permutation);

static
size_t binary_search(const uint8_t *vector, size_t element_size, size_t size,
                     const uint8_t *key);


static size_t boundary_signature(
  const solver *solver,
  const quad_index *quads, uint32_t num_quads,
  uint8_t *signature, uint8_t *valence) {
  signature[0] = num_quads;

  memset(valence, 0, VERTEX_SET_MAX);
  for (uint32_t i = 0; i < num_quads; i++) {
    const vertex_index *quad = solver->quads + 4 * quads[i];
    for (quad_vertex_index j = 0; j < 4; j++)
      valence[quad[j]]++;
  }

  uint8_t max_valence = 3;
  for (vertex_index v = 0; v < solver->max_vertex_id; v++) {
    if (valence[v] > max_valence) max_valence = valence[v];
  }

  uint8_t *histogram = signature + 1;
  for (vertex_index v = 0; v < solver->max_vertex_id; v++) {
    if (valence[v] != 0)
      histogram[valence[v] - 3]++;
  }

  uint8_t *edge_histogram = histogram + (max_valence - 2);
  for (uint32_t i = 0; i < num_quads; i++) {
    const vertex_index *quad = solver->quads + 4 * quads[i];
    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_index v = quad[j], u = quad[(j + 1) % 4];

      if (valence[v] >= valence[u]) {
        size_t offset = ((valence[v] - 2) * (valence[v] - 3)) / 2;
        edge_histogram[offset + valence[u] - 3]++;
      }
    }
  }

  return 1 + (max_valence - 2) +
    ((max_valence - 2) * (max_valence - 3)) / 2 + 1;
}

void component_lower_bound(
  solver *solver,
  quad_index *quads, uint32_t num_quads,
  vertex_index *v_min,
  hex_index *h_min) {
  if (!lower_bounds.loaded || num_quads > lower_bounds.num_quads) {
    *v_min = 0;
    *h_min = (num_quads - 2 + 3) / 4;
    return;
  }

  uint8_t signature[lower_bounds.signature_size];
  memset(signature, 0, sizeof(signature));

  uint8_t valence[VERTEX_SET_MAX] = {0};
  boundary_signature(solver, quads, num_quads, signature, valence);

  size_t i = binary_search(
    lower_bounds.signatures, lower_bounds.signature_size,
    lower_bounds.signatures_num_bytes / lower_bounds.signature_size,
    signature);

  vertex_index permutation[VERTEX_SET_MAX];
  vertex_index inv_permutation[VERTEX_SET_MAX];

  for (size_t offset = i * lower_bounds.signature_size;
       offset < lower_bounds.signatures_num_bytes &&
         memcmp(lower_bounds.signatures + offset, signature, sizeof(signature)) == 0;
       offset += lower_bounds.signature_size, i++) {
    if (compare_boundary(solver, quads,
                         lower_bounds.meshes + 4 * i * lower_bounds.num_quads,
                         num_quads, valence,
                         permutation, inv_permutation)) {
      *v_min = lower_bounds.minima[2 * i + 0];
      *h_min = lower_bounds.minima[2 * i + 1];
      return;
    }
  }

  *v_min = *h_min = 0;
}

static bool check_isomorphism(
  solver *solver, quad_index *quads,
  const vertex_index *ref, uint32_t num_quads,
  quad_index start_quad,
  quad_vertex_index offset, bool flip,
  const uint8_t *valence_a, const uint8_t *valence_b,
  vertex_index *permutation, vertex_index *inv_permutation);

static
bool compare_boundary(
  solver *solver, quad_index *quads,
  const vertex_index *ref, uint32_t num_quads,
  const uint8_t *valence,
  vertex_index *permutation, vertex_index *inv_permutation) {
  quad_index original_quads[num_quads];
  memcpy(original_quads, quads, sizeof(original_quads));

  uint8_t valence_ref[VERTEX_SET_MAX] = {0};

  if (valence) {
    for (size_t i = 0; i < 4 * num_quads; i++)
      valence_ref[ref[i]]++;
  }

  for (quad_index i = 0; i < num_quads; i++) {
    for (quad_vertex_index offset = 0; offset < 4; offset++) {
      for (int flip = 0; flip < 2; flip++) {
        if (check_isomorphism(solver, quads, ref, num_quads,
                              original_quads[i], offset, flip,
                              valence_ref, valence,
                              permutation, inv_permutation)) {
          memcpy(quads, original_quads, sizeof(original_quads));
          return true;
        }
      }
    }
  }

  memcpy(quads, original_quads, sizeof(original_quads));

  return false;
}

static uint32_t find_quad_with_edge(solver *solver,
                                    quad_index *quads, uint32_t num_quads,
                                    vertex_index a, vertex_index b,
                                    vertex_index *opposite);

static bool check_isomorphism(
  solver *solver, quad_index *quads,
  const vertex_index *ref, uint32_t num_quads,
  quad_index start_quad,
  quad_vertex_index offset, bool flip,
  const uint8_t *valence_a, const uint8_t *valence_b,
  vertex_index *permutation, vertex_index *inv_permutation) {
  for (uint32_t i = 0; i < VERTEX_SET_MAX; i++)
    permutation[i] = inv_permutation[i] = VERTEX_SET_MAX;

  uint32_t i_b;
  for (i_b = 0; quads[i_b] != start_quad; i_b++);

  quads[i_b] = quads[0];
  quads[0] = start_quad;

  static const quad_vertex_index flip_perm[] = {0, 3, 2, 1};
  for (quad_vertex_index i = 0; i < 4; i++) {
    vertex_index v_b;
    if (flip)
      v_b = solver->quads[4 * start_quad + flip_perm[(i + offset) % 4]];
    else
      v_b = solver->quads[4 * start_quad + (i + offset) % 4];

    if (valence_b && valence_a[ref[i]] != valence_b[v_b])
      return false;

    if (!extend_permutation(permutation, inv_permutation, ref[i], v_b))
      return false;
  }

  for (uint32_t i = 1; i < num_quads; i++) {
    vertex_index opposite[2];
    uint32_t i_b = find_quad_with_edge(solver, quads + i, num_quads - i,
                                       permutation[ref[4 * i + 0]],
                                       permutation[ref[4 * i + 1]],
                                       opposite);

    quad_index q_b = quads[i + i_b];
    quads[i + i_b] = quads[i];
    quads[i] = q_b;

    for (quad_vertex_index j = 2; j < 4; j++) {
      if (valence_b && valence_a[ref[4 * i + j]] != valence_b[opposite[j - 2]])
        return false;

      if (!extend_permutation(permutation, inv_permutation,
                              ref[4 * i + j], opposite[j - 2]))
        return false;
    }
  }

  return true;
}

static uint32_t find_quad_with_edge(solver *solver,
                                    quad_index *quads, uint32_t num_quads,
                                    vertex_index a, vertex_index b,
                                    vertex_index *opposite) {
  for (uint32_t i = 0; i < num_quads; i++) {
    quad_index q = quads[i];

    quad_vertex_index j;
    for (j = 0; j < 4 && solver->quads[4 * q + j] != a; j++);

    if (j != 4) {
      if (solver->quads[4 * q + (j + 1) % 4] == b) {
        opposite[0] = solver->quads[4 * q + (j + 2) % 4];
        opposite[1] = solver->quads[4 * q + (j + 3) % 4];
        return i;
      }
      else if (solver->quads[4 * q + (j + 3) % 4] == b) {
        opposite[0] = solver->quads[4 * q + (j + 2) % 4];
        opposite[1] = solver->quads[4 * q + (j + 1) % 4];
        return i;
      }
    }
  }

  return num_quads;
}

static inline
size_t lb_log2(unsigned long x);

static
size_t binary_search(const uint8_t *vector, size_t element_size, size_t size,
                     const uint8_t *key) {
  const uint8_t *low = vector;
  size_t full_size = size;

  if (size <= 1) return 0;

  size_t step = size - (1ul<<(lb_log2(size)-1));
  const uint8_t *mid = low + element_size * step;
  if (memcmp(mid, key, element_size) < 0) low += step * element_size;

  size -= step;

  while (size != 0) {
    size /= 2;
    const uint8_t* mid = low + element_size * size;
    if (memcmp(mid, key, element_size) < 0)
      low += element_size * size;
  }

  while ((size_t)(low - vector) < full_size * element_size &&
         memcmp(low, key, element_size) < 0) {
    low += element_size;
  }

  return (low - vector) / element_size;
}

static inline
size_t lb_log2(unsigned long x) {
  if (x <= 1) return 0;
#ifdef __GNUC__
  return (8*sizeof(unsigned long)) - __builtin_clzl(x-1);
#else
  size_t i = 8*sizeof(unsigned long);
  unsigned long bit = 1ul << (i-1);

  while (((x-1) & bit) == 0) {
    bit >>= 1;
    i--;
  }

  return i;
#endif
}

static
size_t maximum_element_size(uint32_t num_boundary_quads, uint32_t num_vertices);

static
bool upper_bounds_look_up_or_insert(
  upper_bounds_table *table,
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *signature,
  uint8_t **location,
  vertex_index *permutation, vertex_index *inv_permutation);

static
bool upper_bounds_look_up_or_insert_no_grow(
  upper_bounds_table *table,
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *signature,
  uint8_t **location,
  vertex_index *permutation, vertex_index *inv_permutation);

static
error_code upper_bounds_grow(upper_bounds_table *table);

error_code init_upper_bounds_table(upper_bounds_table *table,
                                   uint32_t num_boundary_quads,
                                   uint32_t num_vertices) {
  size_t element_size = maximum_element_size(num_boundary_quads, num_vertices);

  table->element_size = element_size;
  table->capacity = (1 << 10);
  table->num_entries = 0;

  table->hash_table = calloc(table->capacity, element_size);

  table->symmetries_size = 0;
  table->symmetries_capacity = 0;
  table->symmetries = NULL;

  return table->hash_table ? SUCCESS : OUT_OF_MEMORY;
}

void release_upper_bounds_table(upper_bounds_table *table) {
  free(table->symmetries);
  free(table->hash_table);
}

typedef struct upper_bounds_entry {
  vertex_index min_v;
  hex_index min_h;
  vertex_index payload[];
} upper_bounds_entry;

bool look_up_upper_bound(solver *solver, quad_index *quads, uint32_t num_quads) {
  if (!solver->ub_data)
    return false;

  uint8_t signature[6 * solver->num_quads + sizeof(uint32_t)];
  uint8_t valence[VERTEX_SET_MAX] = {0};
  memset(signature, 0, sizeof(signature));
  *(uint32_t*)signature = boundary_signature(
    solver, quads, num_quads, signature + sizeof(uint32_t),
    valence);

  size_t real_size = *(uint32_t*)signature + sizeof(uint32_t) +
         4 * num_quads * sizeof(vertex_index) +
         sizeof(upper_bounds_entry) +
         sizeof(size_t) +
    3 * solver->num_vertices * sizeof(vertex_set);
  size_t max_size = ((upper_bounds_table*)solver->ub_data)->element_size;
  if (real_size > max_size)
    printf("element size: %zu vs. %zu\n", max_size, real_size);

  uint8_t *location;

  vertex_index permutation[VERTEX_SET_MAX];
  vertex_index inv_permutation[VERTEX_SET_MAX];

  if (!upper_bounds_look_up_or_insert(
        solver->ub_data,
        solver, quads, num_quads, signature,
        &location, permutation, inv_permutation)) {
    return false;
  }

  upper_bounds_entry *entry = (upper_bounds_entry*)(
    location + sizeof(uint32_t) + *(uint32_t*)location);

  if (entry->min_v > solver->max_vertex_id ||
      entry->min_h > solver->hex_id) {
    if (entry->min_v > solver->max_vertex_id &&
        entry->min_h > solver->hex_id) {
      entry->min_v = solver->max_vertex_id;
      entry->min_h = solver->hex_id;
    }
    return false;
  }

  return true;
}

static
size_t maximum_element_size(uint32_t num_boundary_quads, uint32_t num_vertices) {
  size_t max_q = num_boundary_quads + num_vertices;
  size_t max_signature = 26 * max_q + num_vertices + sizeof(uint32_t);
  return max_signature +
    sizeof(upper_bounds_entry) +
    4 * max_q * sizeof(vertex_index) +
    sizeof(size_t) +
    3 * sizeof(vertex_set) * num_vertices;
}

static
bool upper_bounds_look_up_or_insert(
  upper_bounds_table *table,
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *signature,
  uint8_t **location,
  vertex_index *permutation, vertex_index *inv_permutation) {
  if ((table->num_entries + 1) * 100 >= 80 * table->capacity) {
    if (upper_bounds_grow(table) != SUCCESS)
      return false;
  }

  return upper_bounds_look_up_or_insert_no_grow(
    table, solver, quads, num_quads,
    signature, location,
    permutation, inv_permutation);
}

static
uint32_t hash_entry(const uint8_t *entry);

static
bool compare_with_entry(solver *solver, quad_index *quads, uint32_t num_quads,
                        const uint8_t *signature,
                        const uint8_t *entry,
                        vertex_index *permutation, vertex_index *inv_permutation);

static
void encode_entry(uint8_t *location,
                  solver *solver, quad_index *quads, uint32_t num_quads,
                  const uint8_t *signature);

static
bool upper_bounds_look_up_or_insert_no_grow(
  upper_bounds_table *table,
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *signature,
  uint8_t **location,
  vertex_index *permutation, vertex_index *inv_permutation) {
  size_t bucket = hash_entry(signature) & (table->capacity - 1);
  size_t i = 0;

   while (true) {
    i++;

     if (*(uint32_t*)(table->hash_table + (table->element_size * bucket)) == 0) {
      table->num_entries++;
      *location = table->hash_table + table->element_size * bucket;
      encode_entry(*location, solver, quads, num_quads, signature);
      return false;
    }
    else if (compare_with_entry(
               solver, quads, num_quads, signature,
               table->hash_table + table->element_size * bucket,
               permutation, inv_permutation)) {
      *location = table->hash_table + table->element_size * bucket;
      return true;
    }
    else {
      bucket = (bucket + i) & (table->capacity - 1);
    }
  }
}

static
uint32_t hash_entry(const uint8_t *entry) {
  uint32_t key_size = *(const uint32_t*)entry;
  uint32_t hash;
  MurmurHash3_x86_32(entry + sizeof(uint32_t), key_size, 0, &hash);
  return hash;
}

static
bool compare_with_entry(solver *solver, quad_index *quads, uint32_t num_quads,
                        const uint8_t *signature,
                        const uint8_t *entry,
                        vertex_index *permutation, vertex_index *inv_permutation) {
  uint32_t size_a = *(const uint32_t*)signature;
  uint32_t size_b = *(const uint32_t*)entry;

  if (size_a != size_b) return false;

  if (memcmp(signature + sizeof(uint32_t),
             entry + sizeof(uint32_t), size_a) != 0)
    return false;

  const upper_bounds_entry *value = (const upper_bounds_entry*)(
    entry + sizeof(uint32_t) + size_a);

  return compare_boundary(solver, quads, value->payload, num_quads,
                          NULL, permutation, inv_permutation);
}

static
void encode_entry(uint8_t *location,
                  solver *solver, quad_index *quads, uint32_t num_quads,
                  const uint8_t *signature) {
  memcpy(location, signature, sizeof(uint32_t) + *(const uint32_t*)signature);
  location += sizeof(uint32_t) + *(const uint32_t*)signature;

  upper_bounds_entry *entry = (upper_bounds_entry*)location;

  entry->min_v = solver->max_vertex_id;
  entry->min_h = solver->hex_id;

  size_t offset = 0;

  bool seen[solver->max_quad_id];
  for (quad_index i = 0; i < solver->max_quad_id; i++) seen[i] = false;

  quad_edge queue[num_quads];
  size_t queue_size = 0;

  queue[queue_size++] = make_quad_edge(quads[0], 0);
  seen[quads[0]] = true;

  for (size_t i = 0; i < queue_size; i++) {
    quad_edge edge = queue[i];
    quad_index q = quad_edge_quad(edge);
    quad_vertex_index first = quad_edge_edge(edge);

    for (quad_vertex_index j = 0; j < 4; j++)
      entry->payload[offset++] = solver->quads[4 * q + (first + j) % 4];

    for (quad_edge_index j = 0; j < 4; j++) {
      quad_edge neighbor = solver->quad_adjacency[4 * q + j];
      if (!seen[quad_edge_quad(neighbor)]) {
        seen[quad_edge_quad(neighbor)] = true;
        queue[queue_size++] = neighbor;
      }
    }
  }

  *(size_t*)(entry->payload + offset) = SIZE_MAX;

  vertex_set *sets = (vertex_set*)(entry->payload + offset + sizeof(size_t));

  vertex_set *edges = sets;
  vertex_set *diagonals = edges + solver->num_vertices;
  vertex_set *quad_diagonals = diagonals + solver->num_vertices;

  for (vertex_index v = 0; v < solver->num_vertices; v++) {
    edges[v] = solver->known_neighbors[v];
    vertex_set_subtract(edges + v, &solver->finished_vertices);

    diagonals[v] = solver->known_diagonals[v];
    vertex_set_subtract(diagonals + v, &solver->finished_vertices);

    vertex_set_init_empty(quad_diagonals + v);
  }

  for (uint32_t i = 0; i < solver->max_quad_id; i++) {
    const vertex_index *quad = solver->quads + 4 * i;
    for (quad_vertex_index j = 0; j < 2; j++) {
      if (!vertex_set_contains(&solver->finished_vertices, quad[j]) &&
          !vertex_set_contains(&solver->finished_vertices, quad[j + 2])) {
        vertex_set_insert(quad_diagonals + quad[j], quad[j + 2]);
        vertex_set_insert(quad_diagonals + quad[j + 2], quad[j]);
      }
    }
  }
}

static
void upper_bounds_reinsert(upper_bounds_table *table, const uint8_t *entry);

static
error_code upper_bounds_grow(upper_bounds_table *table) {
  uint8_t *old_table = table->hash_table;
  size_t old_capa = table->capacity;

  uint8_t *new_table = calloc(2 * table->capacity, table->element_size);;
  if (!new_table)
    return OUT_OF_MEMORY;

  table->capacity *= 2;
  table->hash_table = new_table;

  for (size_t i = 0; i < old_capa; i++) {
    if (*(uint32_t*)(old_table + table->element_size * i) == 0)
      continue;
    upper_bounds_reinsert(table, old_table + table->element_size * i);
  }

  free(old_table);
  return SUCCESS;
}

static
void upper_bounds_reinsert(upper_bounds_table *table, const uint8_t *entry) {
  size_t bucket = hash_entry(entry) & (table->capacity - 1);
  size_t i = 0;

  while (true) {
    i++;

    if (*(uint32_t*)(table->hash_table + (table->element_size * bucket)) == 0) {
      memcpy(table->hash_table + table->element_size * bucket,
             entry, table->element_size);
      return;
    }
    else {
      bucket = (bucket + i) & (table->capacity - 1);
    }
  }
}

static
bool yield_best_from_known_solution(
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *valence, size_t sol_id);

void check_known_solutions(solver *solver,
                           quad_index *quads, uint32_t num_quads) {
  if (!known_solutions.loaded)
    return;

  if (num_quads > known_solutions.num_quads)
    return;

  uint8_t signature[known_solutions.signature_size];
  uint8_t valence[VERTEX_SET_MAX] = {0};
  memset(signature, 0, sizeof(signature));
  boundary_signature(solver, quads, num_quads, signature, valence);

  size_t i = binary_search(
    known_solutions.signatures, known_solutions.signature_size,
    known_solutions.signatures_num_bytes / known_solutions.signature_size,
    signature);

  for (size_t offset = i * known_solutions.signature_size;
       offset < known_solutions.signatures_num_bytes &&
         memcmp(known_solutions.signatures + offset, signature, sizeof(signature)) == 0;
       offset += known_solutions.signature_size, i++) {
    if (yield_best_from_known_solution(solver, quads, num_quads, valence, i))
      break;
  }

  return;
}

static
size_t build_from_known_solution(solver *solver, quad_index *quads, uint32_t num_quads,
                                 vertex_index *permutation, bool flip,
                                 size_t sol_id, vertex_index *solution);

static
bool yield_best_from_known_solution(
  solver *solver, quad_index *quads, uint32_t num_quads,
  const uint8_t *valence, size_t sol_id) {
  const vertex_index *ref = known_solutions.meshes +
    4 * sol_id * known_solutions.num_quads;

  vertex_index permutation[VERTEX_SET_MAX];
  vertex_index inv_permutation[VERTEX_SET_MAX];

  quad_index original_quads[num_quads];
  memcpy(original_quads, quads, sizeof(original_quads));

  hex_index num_added = known_solutions.minima[2 * sol_id + 1];
  vertex_index best_solution[8 * (solver->hex_id + num_added + num_quads)];
  size_t min_solution_size = SIZE_MAX;

  bool found_match = false;

  uint8_t valence_ref[VERTEX_SET_MAX] = {0};
  if (valence) {
    for (size_t i = 0; i < 4 * num_quads; i++)
      valence_ref[ref[i]]++;
  }

  for (quad_index i = 0; i < num_quads; i++) {
    for (quad_vertex_index offset = 0; offset < 4; offset++) {
      for (int flip = 0; flip < 2; flip++) {
        if (check_isomorphism(solver, quads, ref, num_quads,
                              original_quads[i], offset, flip,
                              valence_ref, valence,
                              permutation, inv_permutation)) {
          found_match = true;

          memcpy(quads, original_quads, sizeof(original_quads));
          vertex_index solution[8 * (solver->hex_id + num_added + num_quads)];
          size_t solution_size =
            build_from_known_solution(solver, quads, num_quads,
                                      permutation, flip,
                                      sol_id, solution);

          if (solution_size < min_solution_size) {
            min_solution_size = solution_size;
            memcpy(best_solution, solution, 8 * solution_size * sizeof(*solution));
          }
        }
      }
    }
  }

  if (found_match) {
    hex_index *old_hex = solver->hexes;

    solver->hexes = best_solution;
    solver->cb(solver, min_solution_size, solver->user_data);

    solver->hexes = old_hex;
  }

  memcpy(quads, original_quads, sizeof(original_quads));
  return found_match;
}

static
size_t count_hexahedra_to_split(solver *solver,
                                const quad_index *quads, uint32_t num_quads,
                                const vertex_index *solution, hex_index num_hexes,
                                hex_index *to_split, bool *unsplittable);

static
size_t pillow_hexahedra(solver *solver, const hex_index *to_split, size_t num_split,
                        vertex_index *solution, hex_index num_hexes, vertex_index num_vertices);

static
size_t build_prism_layer(solver *solver, quad_index *quads, uint32_t num_quads,
                         vertex_index *solution, hex_index num_hexes, vertex_index num_vertices);

static
size_t build_from_known_solution(solver *solver, quad_index *quads, uint32_t num_quads,
                                 vertex_index *permutation, bool flip,
                                 size_t sol_id,
                                 vertex_index *solution) {
  hex_index num_added = known_solutions.minima[2 * sol_id + 1];

  memcpy(solution, solver->hexes, 8 * solver->hex_id * sizeof(*solution));

  vertex_index max_v = solver->max_vertex_id;

  for (hex_index j = 0; j < num_added; j++) {
    vertex_index hex[8];
    memcpy(hex, known_solutions.hexahedra +
           8 * (sol_id * known_solutions.num_hexes + (num_added - j - 1)),
           sizeof(hex));

    for (hex_vertex_index k = 0; k < 8; k++) {
      if (permutation[hex[k]] == VERTEX_SET_MAX)
        permutation[hex[k]] = max_v++;
      hex[k] = permutation[hex[k]];
    }

    if (!flip) {
      vertex_index tmp;

      tmp = hex[1]; hex[1] = hex[3]; hex[3] = tmp;
      tmp = hex[5]; hex[5] = hex[7]; hex[7] = tmp;
    }

    memcpy(solution + 8 * (solver->hex_id + j), hex, sizeof(hex));
  }

  size_t solution_size = solver->hex_id + num_added;

  hex_index to_split[num_added];
  bool unsplittable;
  size_t incompatibilities = count_hexahedra_to_split(solver,
                                                      quads, num_quads,
                                                      solution, solution_size,
                                                      to_split, &unsplittable);

  if (6 * incompatibilities < num_quads && !unsplittable)
    solution_size = pillow_hexahedra(solver, to_split, incompatibilities,
                                     solution, solution_size, max_v);
  else {
    size_t old_size = solution_size;
    solution_size = build_prism_layer(solver, quads, num_quads,
                                      solution, solution_size, max_v);
    size_t delta_size = solution_size - old_size;

    vertex_index added[8 * num_added];
    memcpy(added, solution + 8 * solver->hex_id, sizeof(added));
    memmove(solution + 8 * solver->hex_id, solution + 8 * (solver->hex_id + num_added),
            8 * delta_size * sizeof(*solution));
    memcpy(solution + 8 * (solver->hex_id + delta_size), added, sizeof(added));
  }

  return solution_size;
}

static
bool find_hex_adjacent_to_facet(solver *solver, const vertex_index *hex, hex_facet_index i,
                                hex_index *out);

static
bool contains_facet_adjacent_to(solver *solver, const vertex_index *hex, hex_index other,
                                hex_facet_index excluded);

static
size_t count_hexahedra_to_split(solver *solver,
                                const quad_index *quads, uint32_t num_quads,
                                const vertex_index *solution, hex_index num_hexes,
                                hex_index *to_split, bool *unsplittable) {
  vertex_set boundary_diagonals[VERTEX_SET_MAX];
  for (size_t i = 0; i < solver->max_vertex_id; i++)
    vertex_set_init_empty(boundary_diagonals + i);

  for (uint32_t i = 0; i < num_quads; i++) {
    const vertex_index *quad = solver->quads + 4 * quads[i];
    for (quad_vertex_index j = 0; j < 4; j++)
      vertex_set_insert(boundary_diagonals + quad[j], quad[(j + 2) % 4]);
  }

  *unsplittable = false;
  size_t count = 0;

  for (hex_index i = solver->hex_id; i < num_hexes; i++) {
    bool must_be_split = false;
    for (hex_facet_index j = 0; j < 6; j++) {
      vertex_index face[4];
      for (quad_vertex_index k = 0; k < 4; k++)
        face[k] = solution[8 * i + HEX_FACET_CONTENTS[j][k]];

      size_t d0 = diagonal_index(solver->num_vertices, face[0], face[2]);
      size_t d1 = diagonal_index(solver->num_vertices, face[1], face[3]);

      if ((face[0] < solver->max_vertex_id &&
           face[2] < solver->max_vertex_id &&
            solver->quad_descriptors[d0] != EMPTY_QUAD_DESCRIPTOR &&
           !vertex_set_contains(boundary_diagonals + face[0], face[2])) ||
          (face[1] < solver->max_vertex_id &&
           face[3] < solver->max_vertex_id &&
           solver->quad_descriptors[d1] != EMPTY_QUAD_DESCRIPTOR &&
           !vertex_set_contains(boundary_diagonals + face[1], face[3]))) {
        *unsplittable = true;
        must_be_split = true;
      }

      hex_index other_side;
      if (!find_hex_adjacent_to_facet(solver, solution + 8 * i, j, &other_side) ||
          other_side == NO_HEXAHEDRON)
        continue;

      if (contains_facet_adjacent_to(solver, solution + 8 * i, other_side, j))
        must_be_split = true;
    }

    for (quad_vertex_index j = 0; j < 4; j++) {
      if (vertex_set_contains(solver->known_diagonals + solution[8 * i + j],
                              solution[8 * i + 4 + (j + 2) % 4]) ||
          vertex_set_contains(solver->known_neighbors + solution[8 * i + j],
                              solution[8 * i + 4 + (j + 2) % 4]))
        must_be_split = true;
    }

    for (hex_vertex_index j = 0; j < 8; j++) {
      for (uint8_t k = 0; k < 3; k++) {
        if (HEX_NEIGHBORS[j][k] < j) {
          vertex_index a = solution[8 * i + j];
          vertex_index b = solution[8 * i + HEX_NEIGHBORS[j][k]];
          if (vertex_set_contains(solver->known_diagonals + a, b)) {
            *unsplittable = true;
            must_be_split = true;
          }
        }
      }
    }

    if (must_be_split)
      to_split[count++] = i;
  }

  return count;
}

static
bool contains_facet_adjacent_to(solver *solver, const vertex_index *hex, hex_index other,
                                hex_facet_index excluded) {
  for (hex_facet_index j = 0; j < 6; j++) {
    if (j == excluded)
      continue;

    hex_index id;
    if (find_hex_adjacent_to_facet(solver, hex, j, &id) &&
        id == other)
      return true;
  }

  return false;
}

static
bool find_hex_adjacent_to_facet(solver *solver, const vertex_index *hex, hex_facet_index i,
                                hex_index *out) {
  vertex_index a = hex[HEX_FACET_CONTENTS[i][0]];
  vertex_index b = hex[HEX_FACET_CONTENTS[i][2]];

  if (a >= solver->num_vertices || b >= solver->num_vertices)
    return false;

  quad_descriptor desc = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, a, b)];

  if (desc == EMPTY_QUAD_DESCRIPTOR)
    return false;

  *out = quad_descriptor_hex_index(desc);

  return true;
}

static
size_t pillow_hexahedra(solver *solver, const hex_index *to_split, size_t num_split,
                        vertex_index *solution, hex_index num_hexes, vertex_index num_vertices) {
  for (size_t i = 0; i < num_split; i++) {
    vertex_index *hex = solution + 8 * to_split[i];

    vertex_index interior[8];
    for (hex_vertex_index j = 0; j < 8; j++) interior[j] = num_vertices++;

    for (hex_facet_index j = 0; j < 6; j++) {
      vertex_index *new_hex = solution + 8 * num_hexes;

      const quad_vertex_index map[] = {0, 3, 2, 1};
      for (quad_vertex_index k = 0; k < 4; k++) {
        new_hex[k] = hex[HEX_ORIENTED_FACET_CONTENTS[j][map[k]]];
        new_hex[4 + k] = interior[HEX_ORIENTED_FACET_CONTENTS[j][map[k]]];
      }

      num_hexes++;
    }

    for (hex_vertex_index j = 0; j < 8; j++) hex[j] = interior[j];
  }

  return num_hexes;
}

static
size_t build_prism_layer(solver *solver, quad_index *quads, uint32_t num_quads,
                         vertex_index *solution, hex_index num_hexes, vertex_index num_vertices) {
  vertex_index permutation[VERTEX_SET_MAX];
  memset(permutation, VERTEX_SET_MAX, sizeof(permutation));

  for (uint32_t i = 0; i < num_quads; i++) {
    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_index v = solver->quads[4 * quads[i] + j];
      if (permutation[v] == VERTEX_SET_MAX)
        permutation[v] = num_vertices++;
    }
  }

  for (size_t i = 8 * solver->hex_id; i < 8 * num_hexes; i++) {
    if (permutation[solution[i]] != VERTEX_SET_MAX)
      solution[i] = permutation[solution[i]];
  }

  for (uint32_t i = 0; i < num_quads; i++, num_hexes++) {
    memcpy(solution + 8 * num_hexes, solver->quads + 4 * quads[i],
           4 * sizeof(*solution));

    for (quad_vertex_index j = 0; j < 4; j++)
      solution[8 * num_hexes + 4 + j] = permutation[solution[8 * num_hexes + j]];
  }

  return num_hexes;
}
